import numpy as np
import scipy as sc
import matplotlib.pyplot as mp
import interpol as it
import load_foil as lf
from scipy import integrate

#---- Methods of approched integrations -----#


def general_method(type_of_method, f, n, a, b):
    "This function calls the adapted method"
    return type_of_method(f, n , a, b)


def left_rectangle_method(f, n, a, b):
    "This function codes the integrating method by left rectangles"
    h = (b - a) / n
    tot = 0

    for i in np.arange(1, n+1, 1):
        tot = tot + f(a + i*h)

    return h*tot


def right_rectangle_method(f, n, a, b):
    "This function codes the integrating method by right rectangles"
    h = (b - a) / n
    tot = 0
    
    for i in range(n):
        tot = tot + f(a + i*h)

    return h*tot


def left_trapezium_method(f, n, a, b):
    "This function codes the integrating method by left trapeziums"
    h = (b - a) / n
    tot = 0

    tot = left_rectangle_method(f, n, a, b) - h/2 * (f(a) - f(b))
    return tot


def right_trapezium_method(f, n, a, b):
    "This function codes the integrating method by right trapeziums"
    h = (b - a) / n
    tot = 0

    tot = right_rectangle_method(f, n, a, b) - h/2 * (f(b) - f(a))
    return tot

def middle_dot_method(f, n, a, b):
    "This function codes the integrating method by middle dots"
    h = (b - a) / n
    tot = 0
    
    for i in range(n):
        tot = tot + f(a + i*h + h/2)

    return h*tot

def simpson_method(f, x, y):
    "This function codes the simpson method"
    x1=x
    y1=f(x1)
    x2=(x+y)/2.0
    y2=f(x2)
    x3=y
    y3=f(x3)
    
    a2 = ((y3-y2)/(x3-x2) - (y2-y1)/(x2-x1))/(x3-x1) 	# calcul du coefficient de x^2
    a1 = ((y2-y1) - a2*(x2*x2 - x1*x1))/(x2-x1)		# calcul du coefficient de x
    a0 = y1 - a2*x1**2 - a1*x1				# calcul du coefficient constant
    
    u = (a2/3.0)*(x1**3) + (a1/2.0)*(x1**2) + a0*x1		# calcul de la primitive en a
    v = (a2/3.0)*(x3**3) + (a1/2.0)*(x3**2) + a0*x3		# calcul de la primitive en b
    
    return (v-u)


def call_simpson(f, N, a, b):
    "This function calls the simpson method for each interval"
    h = (b - a) / N
    tot = 0
    for i in range(n):
        tot = tot + simpson_method(f, a + i*h, a+(i+1)*h)
    return tot


# A partir des differents tests de convergence effectues dans converge_tests.py, il s'avere que malgre notre implementation
# ce n'est pas la methode de Simpson mais celle du point milieu qui converge le plus vite.
# Nous avons donc choisi dans la suite d'utiliser celle-ci.

# ----- Functions to compute the derivate ----- #

def derivate_function(ex,ey,y2,j):
    return lambda x : (((ey[j+1]-ey[j])/(ex[j+1]-ex[j]))-((3*((ex[j+1]-x)/(ex[j+1]-ex[j]))**2-1)/6)*(ex[j+1]-ex[j])*y2[j]+((3*((x - ex[j])/(ex[j+1]-ex[j]))**2-1)/6)*(ex[j+1]-ex[j])*y2[j+1]) 

    
def derivate(y2,n,ex,ey):
    f1 =  []
    for j in range(n-1):
       
       f1 += [derivate_function(ex,ey,y2,j)]
     
    return f1

print "test derivate function :"
f = derivate(np.array([5,6,4]),3,np.array([2,3,1]),np.array([1.5,0.4,1]))
print f[0](5)
print f[1](5)


print "general_method(middle_dot_method,lambda x: 2*x,6,0,2) : "
print general_method(middle_dot_method,lambda x: 2*x,6,0.,2.)
print sc.integrate.fixed_quad(lambda x: 2*x,0.,2.)



y2 = it.spline(it.ex,it.ey,it.ne,it.yp1,it.ypn)
n = it.ne



#----- Function to integer ------#

def integrale(y2,n,ex,ey,method):
    length = 0.
    #test = 0.
    d = derivate(y2,n,ex,ey)
    
    for i in np.arange(0,np.shape(d)[0],1):
        der = d[i]
        F = lambda x: np.sqrt(1+der(x)**2)

        length += general_method(method,F,n,ex[i],ex[i+1])
        #test += sc.integrate.fixed_quad(F,ex[i],ex[i+1])[0]
        #print length
        #print test
        
    return length

print integrale(y2,n,it.ex,it.ey,middle_dot_method)
